Multiband bosons in optical lattices 
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We study a gas of repulsively interacting bosons in an optical lattice and explore the physics 
beyond the lowest band Hubbard model. Utilizing a generalized Gutzwiller ansatz, we find how 
the lowest band physics is modified by the inclusion of the first excited bands. In contrast to the 
prediction of the lowest band Bose-Hubbard model, a reentrant behavior of superfluidity is envisaged 
as well as decreasing width of the Mott lobes at strong coupling. 
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I. INTRODUCTION 

Experiments with cold atoms in optical lattices have 
made dramatic progress in recent past [J,!!. Optical lat- 
tices, low densities, and low temperatures render fantas- 
tic degree of control and have enabled detailed studies of 
strongly correlated quantum svstems. For example, the 
Mott-superfluid transition 0, has been successfully ob- 
served in optical lattices. This transition can occur even 
at T = and is due to the competition between kinetic 
energy and repulsive on-site interactions between lowest 
band bosons. For large interactions, the energy is min- 
imized in an incompressible state with fixed atom num- 
bers while for weaker interactions kinetic energy favors 
atomic tunneling and drives the system into a superfluid. 

Using Feshbach resonances [5] , or by deepening the lat- 
tice, one can increase the strength of the atom-atom in- 
teraction compared with the kinetic energy of the atoms 
in the lattice. These stronger interactions can cause 
considerable excited band populations [6] and take the 
system beyond the region where the usual lowest band 
description is adequate for a quantitatively accurate 
description 01 • Substantial excited band populations 
caused by interactions have indeed been experimentally 
observed 0]. In another interesting experiment, Winkler 
et al. used Feshbach resonance to increase the interac- 
tion energy so much that the energy of two atom states 
lied in the bandgap. In this way they demonstrated the 
existence of repulsively bound pairs, which were unable 
to dissociate due to the absence of available final atomic 
states. However, these experiments did not probe the re- 
gion where the energy would have been even higher and 
dissociation to excited bands would have become possi- 
ble. 

Another way for the excited bands to play a relevant 
role is to directly couple atoms from the lowest band 
to the excited bands. This was experimentally demon- 
strated recently by Miiller et al. [1^ by coupling atoms 
from the lowest band Mott insulator into the first ex- 
cited p-band of the lattice via Raman transitions between 
bands. They found lifetimes of p-band atoms which were 



considerably longer than the tunneling time-scale in the 
lattice and were also able to explore how coherence on 
the excited band was established. This experiment paves 
the way to explore the equilibrium physics of the purely 
p-band bosons |11] and outlines a pos sible route to real- 
ize supersolids [12] or novel phases [l3',^l^ on the excited 
bands of an optical lattice. Alternatively, higher bands 
can be populated by fermions if the filling factor is larger 
than one [ISj]. In this case the Pauli exclusion principle 
ensures that the Fermions that cannot populate the low- 
est band, must occupy the excited bands [16] . When the 
optical lattice geometry is varied from the usual cubic 
structure, exciting analogs also appear between excited 
orbital physics and graphene [17|]. 



In this paper we explore the role of the first excited 
bands on the phase diagrams of the Bose-Hubbard model. 
We will show that the Mott lobes can become strongly 
deformed due to combined effects of interactions and ex- 
cited bands. In some parameter regions, reentrant behav- 
ior of the insulator superfluid transitions may arise and 
in the strongly interacting regime, complex phase dia- 
grams with the possibility of excited band superfluidity 
appears. 
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The paper is organized as follows. In Sec.[TTlwe discuss 
the physical system and derive the Hamiltonian of the 
system. Here we also explain the variational approach we 
use to solve for the ground state configuration. In Sec. lIIII 
we proceed by solving for the multiband phase diagrams 
in different dimensions and discuss the salient features of 
how excited bands change the lowest band physics. We 
also discuss briefly how our results would relate to results 
in a trap and what kind of changes excited band atoms 
will cause in the trapped system compared with the usual 
lowest band results. We conclude with some discussion 
in SecIVl 
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II. THEORY AND AND GENERALIZED 
GUTZWILLER ANSATZ 

The microscopic Hamiltonian for the dilute Bose gas 
at low temperatures in a trap is given by 



and 



microscopic 



2m 



V{v) 



^t/;t (r)t/;t (r)^(r)^(r) - (r)^(r) , 

(1) 

where fi is the chemical potential, m the atomic mass, g is 
the interatomic interaction strength, and ?/^(r) and ip^r) 
are the bosonic annihilation and creation operators, while 
V{v) is the external trapping potential which in this work 
is taken to be the lattice potential 



V(r) =Vl sin^ 



I d /' 



(2) 



where d is the lattice spacing and Vl the lattice depth. 
For a deep lattice it is reasonable to expand field oper- 
ators in terms of the localized Wannier states. Here we 
go beyond the usual lowest band Hubbard model by also 
including the first excited states (p-band). In a three di- 
mensional lattice this implies an expansion of the field 
operators 



(3) 



where i = (ix^iy^iz) labels the lattice site and a G 
{0, X, z} is the fiavor index. The bosonic operators 
annihilate a boson of flavor a from the site i. We com- 
pute the Wannier functions and bandgaps from the ideal 
gas band structure. 

Substituting the operator expansions into Eq. ^ and 
ignoring all but the leading order on-site interactions and 
nearest neighbor tunneling processes we derive our fun- 
damental Hamiltonian 



H = Ho + Hr, 



where the ideal part is given by 



Hfd, 



(4) 



1 a, a <i,j>a 

(5) 

Here ABG,a are the bandgaps and Xl^i indicates the 
nearest neighbor sum over the neighbors in the direction 
a G {x^y^z}. Since the Bloch functions diagonalize the 
one-body hamiltonian, there are no interband hopping 
terms in the Wannier representation considered here [18j. 
The terms originating from the interatomic interactions 
are given by 

i cr i ^a^(j' 

(6) 



(J(j' ^<J^<J' 



where Hfd contains terms that describe flavor chang- 
ing collisions and collisions transferring atoms between 
bands. This term has formal similarity with terms re- 
sponsible for spin-dynamics in spinor condensates [l^, 
zOj. However, the strength of these terms is comparable 
to other interaction terms as opposed to spinor conden- 
sates where it is usually small, being proportional to the 
difference between singlet and triplet scattering lengths 
(for spin-1 spinor condensate). If the lattice site were pre- 
cisely harmonic there could also be processes where two 
atoms on the p-band collide and scatter into an atom 
on the lowest band and an atom on the d-band. In a 
real lattice, such processes are off-resonant and they are 
also ignored in our model restricted to just the lowest 
bands. For the lattice depth considered in this paper, 
the anharmonicity is indeed evident, and only momentum 
states close to the Brillouin edges are non far off-resonant. 
In addition, the population of p-band bosons is in gen- 
eral small and the corresponding scattering amplitudes 
have been verified to be typically one order of magnitude 
smaller than for non-flavor changing collisions. 

Further, terms describing scattering between atoms in 
neighboring sites have been left out as well. For our 
choice of lattice depth, we have checked that the mag- 
nitude of these particle assisted tunneling terms are 1-2 
% of the 5-band on-site scattering amplitude, and hence 
such an approximation is justified. It should be kept in 
mind though, that for other system parameters there are 
circumstances when nearest neighbor interactions [12] or 
particle assisted tunneling processes [21] might give rise 
to new physics, but these regimes are not considered here. 

The various coupling strengths in the lattice model are 
related to g through 



'I 



(8) 



and the tunneling coefficients are given by 



drwa{r) 



2m 



Wa{r-dea), (9) 



where is the unit vector in the direction a. When 
the lattice is symmetric, the tunneling strength on the 
lowest band is independent of direction. However, the 
directional dependence of the tunneling strength must 
be kept for the p-band atoms, since the overlap integrals 
are very different depending on whether one is integrat- 
ing along the node of the Wannier function or orthogonal 
to it. This can have important consequences for the char- 
acteristics in these systems [ll|, |22| 

The above formulation was derived for the three- 
dimensional situation with 4 fiavors. In this paper we 
will, however, also consider the corresponding special 
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cases of one- and two-dimensional systems with 2 and 
3 flavors respectively. The number of relevant flavors can 
be reduced in asymmetric lattices, where the bandgap 
in the direction of deep lattice potential becomes higher 
than in other directions. 

We approach the physics of the above theory using a 
Gutzwiller ansatz for the many body wave-function gen- 
eralized to multiple flavors [23]. Our ansatz is given by 

l^> = nE/nV)i, (10) 
i n 

where /n ^ are the variational parameters to be deter- 
mined by minimizing the energy. In our Gutzwiller ap- 
proach, the on-site states are expanded in terms of the 
Fock states |n) = |no, n^^, n^, n^) of the multiple flavor 
system. For numerical reasons the sum over n must be 
cut-off and in our computations we include all the states 
with a total excitation Ua- less than 8. The num- 
ber of variational parameters in a site increase propor- 
tional to n^^^^jj where Nf is the number of flavors. For 
this reason, the problem is substantially heavier in terms 
of computing power than implementing the Gutzwiller 
ansatz for a single flavor problem. 

For the lowest band Bose-Hubbard model it is known 
that the Gutzwiller method provides a reasonably accu- 
rate description of the transitions from Mott-insulators 
in three- and two-dimensions. In a more strongly corre- 
lated one-dimensional system it is reasonable only qual- 
itatively with considerable quantitative differences from 
exact computations. This gives us confidence that the 
above ansatz is able to capture most of the interesting 
physics expected to occur in our model. 

We minimize the energy functional at T = using a 
conjugate gradient method. Since the signs of the tun- 
neling strengths can vary and flavor changing collisions 
are sensitive to phase factors in the wave-function am- 
plitudes, we have observed that one has to be careful in 
ensuring a convergence to the real lowest energy state as 
opposed to some local minima. Among other things, neg- 
ative tunneling strength in the kinetic energy can favor tt- 
phase modulation of the condensate order parameter and 
if this modulation is not present in the initial guess input, 
the numerics might not converge to the real minimum. 
Flavor dynamics is also sensitive to the phase factors be- 
tween flavors and it is easy to converge to solutions with 
incorrect relative phase factors between different flavor 
order parameters. 

III. RESULTS 

A. One- dimensional lattice with two flavors 

In a one-dimensional lattice our model reduces to the 
lowest band and one excited p-band. We summarize our 
results in Figs. [1] and O Throughout, the figures dis- 
play the order parameters (V^a), on-site total and lowest 



band number fluctuations An|. and Atiq, and on-site to- 
tal number and flavor number of atoms ut and n^- These 
quantities are presented as functions of the scaled hop- 
ping parameter zt/Uoo (where z is the number of nearest 
neighbours and t the s-band tunneling coefficient) and 
the scaled chemical potential /i//7oo- The lattice depth 
for ^^Rb atoms is fixed at Vl = where En is the 

recoil energy. 

When there is only one atom per site, this atom resides 
on the lowest band and, for relatively weak coupling, 
also for higher atom numbers the result is essentially the 
same as for the lowest band Bose-Hubbard model. As ex- 
pected, one finds a superfluid region (on the lowest band) 
at weak coupling and Mott insulating lobes of different 
integer atom number at somewhat stronger coupling. 

However, at stronger coupling the population of the 
excited band increases smoothly when the on-site atom 
number tit = 2. The density difference between flavors 
starts to fluctuate and the Mott-lobes become distorted 
as they curve toward smaller chemical potentials, giving 
rise to re-entrant behavior not present in the lowest band 
Hubbard model within the Gutzwiller approximation. As 
one moves towards weaker coupling in this region, one 
can start from an insulator, move into a superfluid, then 
back into an insulator and finally re-enter the superfluid 
region. Note, however, that re-entrant behavior could 
exist in one-dimensional models which move beyond the 
Gutzwiller approximation [23]. In that case re-entrant 
behavior is due to correlation effects between sites and is 
not caused by the presence of excited bands. 

At even stronger couplings, the population on the p- 
band becomes more dominant and influences the state 
of the system more strongly. This is depicted in Fig. [21 
which zooms on the regime of small t/UoQ. In the mul- 
tiflavor problem Mott regions can change also qualita- 
tively. In particular, one can often have insulating regions 
where the total on-site atom number is non-fluctuating, 
but where the densities of individual flavors, or their den- 
sity differences do fluctuate. This is possible even in the 
limit of infinitely deep lattice due to flavor changing on- 
site collisions. Physics on the excited band is quite sen- 
sitive to the details of the model, such as the relative 
magnitude of the bandgap and the lowest band band- 
width. With the parameter values than those used here, 
there is a possibility of narrow regions of excited band 
superfluidity appearing in the phase diagram. 

Interestingly there is a peculiar transition around 
fi/Uoo ^ 1.5 at strong coupling. As one goes across this 
region in the direction of weaker coupling one moves from 
a Mott-insulator with ut = 3 predominantly on the ex- 
cited band, into a narrow excited band superfluid region, 
and finally into a Mott-insulator with tit = 3 mainly on 
the lowest band. Across this transition the flavor densi- 
ties (ncr) change abruptly. The difference in the Mott- 
insulators around this transition is that at strong cou- 
pling the state is predominantly of type a|0,3) + /3|2,l) 
while at weaker coupling it is of type a|3,0) + /3|1,2). 
Close to the transition (in the narrow SF region) there 



4 



(a) : KVn)! 



(b) : A n; 




FIG. 1: (Color online) Properties of the one- dimensional two 
flavor Bose-Hubbard model as a function of chemical poten- 
tial and the inverse interaction strength 2t/[/oo, where 2 is 
the number of nearest neighbors. For concreteness the pa- 
rameters were computed for a lattice of depth Vl = 15 Er. 
(a) Lowest band condensate order parameter, (b) total atom 
number fluctuation, (c) fluctuation of the lowest band atom 
number, and (d) the total on-site atom number. The horizon- 
tal yellow line indicates the position t/ Abg and warm (cold) 
colors imply high (low) values. 




FIG. 2: (Color online) Properties of the one-dimensional two 
flavor Bose-Hubbard model as a function of chemical potential 
and the inverse interaction strength 2t/Uoo at strong coupling. 
(a)-(b) Lowest and excited band condensate order parameters, 
(c) total atom number fluctuation, (d)-(e) lowest and excited 
band atom numbers, and (f) the total on-site atom number. 
The horizontal yellow line indicates the position t/ABG x 
10^ and warm (cold) colors imply high (low) values. Lattice 
parameters are the same as in the Fig. [T] 



are more non-zero amplitudes in the wave-function ex- 
pansion. 

In a mean-field Gross-Pitaevskii picture one could 
think of this transition as somewhat analogous to a tran- 
sition in a system of two-component Bose-Einstein con- 
densates. When there is a strong repulsion between com- 



ponents, it is more favorable to put more atoms into the 
state with a maximum density near the edges and a van- 
ishing density in the center. At weaker, coupling it is 
more favorable to have more atoms in the "core" state. 
Related transitions are predicted by the one-dimensional 
multi-orbital mean- field computations [25[. 



B. Two-dimensional lattice with three flavors 

In a two-dimensional symmetric lattice there are two- 
degenerate p-bands which must be taken into account. 
Our results for the two dimensional system are shown in 
Figs. [3] and HJ where [H gives a close-up in the regime of 
small t/UoQ. As can be seen from the figures, the Mott 
lobes are deformed by the excited bands in a similar way 
as in a one-dimensional system with fewer flavors. At 
stronger coupling the lowest band atom number can start 
to fluctuate when the on-site atom number is greater than 
one and this is caused by the possibility of the excited 
band carrying part of the site occupation. 

In the regime of strong coupling, the expected phase di- 
agram is quite complex with the possibility of transitions 
between lowest band and excited band superfluids. Also, 
around /i/t/oo ^ 1-2 at strong coupling there is a tran- 
sition (as one changes t/Uoo) between a Mott insulator 
with nT = 3 to a Mott insulator with ut = 2. Here the 
population of the lowest band is reduced abruptly as one 
moves towards strong coupling. Around fi/Uoo ^ 1-5, a 
similar transition is encountered between insulators with 
nT = 3 and tit = 4 with a narrow superfluid region on 
the excited bands in between the insulating phases. The 
magnitude of the excited band condensate order param- 
eters were found to be equal. 



C. Three-dimensional lattice with four flavors 

In a three-dimensional lattice where the lattice depth 
is equal in all directions, the degeneracy of the p-band is 
increased to 3 and the total number of relevant flavors is 
4. As the dimensionality increases the computations be- 
come substantially more time-consuming partly due to 
larger state space, but also because of more stringent 
convergence properties to the global energy minima for 
large couplings. However, as we demonstrate in Fig. [5] 
the general structure of the phase diagrams remains for 
couplings weaker or of same order of magnitude than 
the bandgaps. The re-entrant behavior of the superfluid- 
insulator transition is therefore independent of dimen- 
sionality. 

For three-dimensional p-band superfluids there is the 
interesting possibility of frustration of the relative phases 
of the condensate order parameters In this paper 

we have assumed that the system is homogeneous, but 
it would be interesting, although computationally more 
challenging, to investigate effects of possible phase frus- 
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FIG. 3: (Color online) Properties of the two-dimensional three 
flavor Bose-Hubbard model as a function of chemical potential 
and the inverse interaction strength 4t/[/oo- The parameters 
are computed for a lattice of depth Vl = VoEr. In (a) we 
show a contour plot of the lowest band condensate order pa- 
rameter. In (b) and (c) we show a contour plot of the total 
number fluctuation and the lowest band number fluctuation, 
respectively. Finally, (d) and (e) show a surface plot of the 
individual flavor densities and (f) the total density. The hor- 
izontal yellow line indicates the position t/A^G and warm 
(cold) colors imply high (low) values. 



FIG. 5: (Color online) Properties of the three-dimensional 
four flavor Bose-Hubbard model as a function of chemical 
potential and the inverse interaction strength 6t/[/oo- The 
parameters are computed for a lattice of depth Vl = \^Er. 
In (a)-(c) we show the lowest band condensate order param- 
eter, total number fluctuation and the number fluctuation of 
the lowest band atoms. Figures (d)-(f) show the lowest band 
atom number, excited band atom number and the total atom 
number respectively. The horizontal line indicates the posi- 
tion i I IS.bg and warm (cold) colors imply high (low) values. 



tration in an inhomogeneous system at very strong cou- 
pling. 
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FIG. 4: (Color online) Properties of the two-dimensional three 
flavor Bose-Hubbard model as a function of chemical potential 
and the inverse interaction strength 4t/?7oo at strong coupling. 
In (a)-(c) we show the lowest and excited band order param- 
eters as well as the total number fluctuation. In (d)-(f) we 
show the lowest and excited band populations together with 
the total on-site atom number. The horizontal yellow line 
indicates the position 10^ x t/A^G and warm (cold) colors 
imply high (low) values. The lattice parameters are the same 
as in Fig. O 



D. Trapped physics 

In the current contribution we have focused on the 
physics of homogeneous systems, but some qualitative 
conclusions about the trapped case are easy enough to 
draw in the spirit of a local density approximation in 
which the chemical potential is replaced by the local 
chemical potential /ie// = M ~ "^(r), where F(r) is the 
trapping potential. For the usual harmonic trap this ef- 
fective chemical potential is maximized in the center of 
the trap and is reduced as one moves to the edge of the 
cloud. 

With respect to our phase diagrams this implies that 
as one moves from the center to the edge of the cloud, 
one moves from right to left on our x-axis /i//7o for some 
particular value of ztjU^^ where z is number of near- 
est neighbors. In the lowest band Bose-Hubbard model 
this suggest either a superfluid everywhere or a wedding 
cake phase structure (at stronger coupling), with Mott- 
insulating plateaus with different atom numbers layered 
in between superfluid regions. 

Inspecting our phase diagrams one can see that this 
picture changes with the inclusion of the excited bands 
in a sense that one can see gradually increasing excited 
band population as the coupling grows. As one increases 
the coupling strength the width of the insulating regions 
first grows, just like in the lowest band Bose-Hubbard 
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model, but eventually starts to decrease, in contrast to 
the predictions of the lowest band Hubbard model. At 
even stronger couplings, a superfluidity on the excited 
band might be established before moving to a fully insu- 
lating phase as coupling approaches infinity. 



IV. DISCUSSION AND CONCLUDING 
REMARKS 

Having found effects due to excited bands in the con- 
text of the bosonic Hubbard model, from an experimental 
point of view, one would naturally wish to measure the 
physical properties of the system in a way which is sen- 
sitive to the band index. In the pioneering experiments 
by Greiner et al. [4] the Mott insulating regime was ob- 
served by vanishing interference peaks after time of flight 
expansion for deep enough lattices. In our model part of 
the systems population lies on the excited bands which 
are typically not phase coherent. After a time of flight 
such contributions from the excited bands would show up 
as incoherent background whose weight depends on the 
fraction of atoms on the excited bands. Furthermore, the 
detailed structure of the interference peaks depends on 
the initial Wannier state prior to expansion and the pres- 
ence of nodes on the Wannier function could, in principle, 
be observed by the careful observation of the structure of 
the interference peaks. However, this might be compli- 
cated in practice especially in the mult i- dimensional sys- 
tems. Addressing the role of flavor dynamics and number 
fluctuation in different flavors requires probes which ad- 
dress different flavors individually. This has been demon- 
strated recently [1^] by coupling atoms with Raman tran- 
sitions between the lowest band and some excited band. 
In Ref. [l^ the selectivity was achieved by breaking the 
degeneracy of p-band atoms by lattice anisotropy, but 
alternatively one could engineer a vanishing transition 
matrix elements for undesired transitions. Interestingly, 
Mueller et al. [10] were also able to observe time-evolution 
of the Bose gas caused by flavor dynamics. 

An additional practical challenge is that most clear de- 
viations from the simple lowest band model occur in the 
regime where the tunneling strength is small compared 
to the interaction strengths. While the strong interac- 
tions will rapidly reach equilibrium in each individual 
site, the establishment of global phase coherence between 
sites would require a longer time since its time-scale is re- 
lated to the tunneling energy. In order to have time to 
establish global phase coherence the onsite loss processes 
should remain sufficiently weak. However, if the interest 
is primarily on on-site properties such as flavor dynamics 
and number fluctuations between flavors, then the re- 
quirements on the strengths of loss processes are not as 
stringent since the relevant time-scale is set by the onsite 
interaction energies. 

Here we solved the multiband Bose-Hubbard model 
utilizing a generalized Gutzwiller ansatz. Having out- 
lined how the couplings to the excited band can change 



the physics from the lowest band Bose-Hubbard model, 
we now proceed to briefly discuss limitations of our for- 
malism. We cut-off our theory to include only the lowest 
band and the first excited band(s) of the lattice. This is 
expected to be a good approximation when the coupling 
is clearly smaller than the bandgaps. However, some of 
our results are computed in a more strongly coupled re- 
gion and in such case one can expect a non-negligible 
level populations on even higher bands [6] and because 
of this, part of our results might change quantitatively 
although not necessarily qualitatively, when more bands 
are included in the theory. However, it is unclear how 
the details of the phase diagrams in the very strongly 
coupled region change if more bands are included in the 
theory. Our main interest has been to elucidate how the 
lowest band results are changed by the inclusion of more 
bands and we have found how the Mott lobes change 
as the couplings become comparable to bandgaps. This 
observation is quite general and not strongly model de- 
pendent. However, details of the excited band physics 
are expected to be more sensitive to the parameters and 
models used. 

With strong interactions and more than one atom per 
site, the localized Wannier functions can be dressed by 
the interactions. Repulsive interactions would spread the 
on-site wave function and therefore lower the on-site in- 
teraction energy and increase the tunneling strength from 
those values used in this paper by computing the ideal 
gas Wannier functions. Such corrections can give rise 
to quantitative corrections [26^ although such corrections 
are not expected to be extremely large for lattices with 
relatively few atoms in each site. Here we have not spec- 
ified how the on-site interaction Uqq is to be related to, 
for example, the three dimensional scattering length a. 
We left this relation unspecified partly because such re- 
lation can depend strongly on dimensionality and partly 
because in the region of very strong coupling the two- 
body scattering physics might change due to confinement 
effects. In particular, if one were to estimate the three- 
dimensional on-site interaction using the usual relation 
^00 = 9 S dY\wQ{Y)\^ together with the harmonic approx- 
imation one can show that Uoo/Abg = \/2/7r (a/a) , 
where a is the size of the Wannier wave-function. This 
shows that interaction becomes comparable to bandgaps 
in the regime where scattering length is comparable to 
the size of the localized wave function and consequently 
the two-body scattering problem might have to be consid- 
ered in the presence of the external potential [27]. If one 
does the similar computation in a two-dimensional pan- 
cake system of thickness a^, one finds Uoo/^bg ^ a/az 
and the coupling becomes comparable to the bandgap 
only in the regime of two-dimensional scattering the- 
ory 0. 

In this paper we have considered the usual single color 
optical lattice and we found that rather strong inter- 
actions are required to induce substantial changes in 
the phase diagram predicted by the lowest band Bose- 
Hubbard model. However, it should be realized that 
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many of the phenomena we discuss here are expected 
to play a role also in a more complicated setting of two- 
color superlattices [2^, [sOI . In such systems one can have 
a lattice of " deep" sites which in turn are effectively dou- 
ble well potentials. In that case the on-site wave function 
on the lowest band is the symmetric wave function which 
peaks in both sites of the deep site. The first excited wave 
function is the antisymmetric and is separated from the 
lowest state by an energy proportional to the tunneling 
strength in the double well system. This bandgap is tun- 
able and can be used to bring the bandgap down even 
when the lattice of deep sites remains well described by 
the tight binding Hubbard model. 

In such system one can expect more dramatic changes 
in the phase diagram at lower interaction strengths. Fur- 
thermore, in superlattices one could engineer band struc- 
tures (in a one-dimensional system) with only two bands 
close in energy, while the next band is substantially sepa- 



rated from the two lowest bands. Hence, under such cir- 
cumstances the theory can be safely cut-off to just include 
the two lowest bands when interactions and bandgap be- 
tween the lowest and the first excited band are smaller 
than the bandgap between first and second excited bands. 
However, the physics might not always be entirely sim- 
ilar to that studied in this paper since density assisted 
tunneling processes and nearest neighbor interactions can 
sometimes play a larger role in superlattices @, ^ . 
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